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Abstract. We show that in a certain, angle-averaged squeezed limit, the IV-point function of 
matter is related to the response of the matter power spectrum to a long-wavelength density 
perturbation, P~ 1 d n P(k\dL)/dd 1 l\s L =o, with n = N — 2. By performing N-body simulations 
with a homogeneous overdensity superimposed on a flat Friedmann-Robertson-Lemaitre- 
Walker (FRLW) universe using the separate universe approach, we obtain measurements of 
the nonlinear matter power spectrum response up to n = 3, which is equivalent to measuring 
the fully nonlinear matter 3— to 5—point function in this squeezed limit. The sub-percent to 
few percent accuracy of those measurements is unprecedented. We then test the hypothesis 
that nonlinear IV-point functions at a given time are a function of the linear power spectrum 
at that time, which is predicted by standard perturbation theory (SPT) and its variants that 
are based on the ideal pressureless fluid equations. Specifically, we compare the responses 
computed from the separate universe simulations and simulations with a rescaled initial 
(linear) power spectrum amplitude. We find discrepancies of 10% at k ~ 0.2 — 0.5/iMpc _1 
for 5— to 3—point functions at 2 = 0. The discrepancy occurs at higher wavenumbers at 
z = 2. Thus, SPT and its variants, carried out to arbitrarily high order, are guaranteed to 
fail to describe matter IV-point functions (N > 2) around that scale. 
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1 Introduction 

Mode-coupling plays a fundamental role in cosmology. Long and short wavelength modes are 
coupled by nonlinear gravitational evolution of matter density fluctuations and the formation 
of dark matter halos and galaxies, as well as via the physics of inflation, i.e., primordial non- 
Gaussianity. The traditional way of characterizing the mode coupling is to study IV-point 
correlation functions of matter, halo, and galaxy density fields, as well as of fluctuations 
of the cosmic microwave background, where the mode coupling enters at N > 2. The iV- 
point functions of matter density fields in particular contain a rich amount of information 
on gravity and the expansion history of the Universe. At N >2, gravity produces specific 
non-Gaussian signatures that can be used to test the physics of structure formation in the 
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Figure 1: Sketch of the squeezed limit configuration of matter A-point functions considered 
in this paper, ki, • • • , k n denote the long-wavelength modes which are spherically averaged 
in Eq. (1.1), while k, k' denote the small-scale modes which are allowed to be fully nonlinear. 


Universe, and must be accounted for when attempting to connect large-scale structure with 
the statistics of the initial conditions in order to search for primordial non-Gaussianity. 

A major challenge is that, beyond the power spectrum, the calculation of (connected) 
nonlinear Appoint functions becomes increasingly difficult. Consider the gold standard for 
predicting the statistics of matter density, N-body simulations. First, a large number of 
modes is necessary to measure the correlation functions with sufficient precision since the 
cosmic variance noise is significant. This demands large computational resources. Second, 
higher A-point functions are more sensitive to transients [1, 2] from the finite starting redshift 
and to mass resolution effects. Third, estimators for higher A-point functions become more 
computationally intensive and difficult to handle. On the theoretical side, the predictions for 
the bispectrum and higher A-point functions likewise become more cumbersome, with terms 
at any given order in perturbation theory rapidly proliferating with A. An analogous effect 
happens in the halo model, with 1- to A-halo terms having to be calculated for the A-point 
function. 

This provides the motivation to study a certain limit considered in this paper where the 
nonlinear A-point functions become simpler and physically more transparent. The limit we 
consider is a specific case of the so-called “squeezed limit”, where there is a hierarchy between 
two large wavenumbers k, k' and A — 2 small wavenumbers ki, ■ ■ • k n . The configuration 
corresponding to this limit is illustrated in figure 1. 

The squeezed limit of dark matter A-point functions has recently been the subject 
of a large body of work in the context of the so-called “consistency relations” [3-15]. The 
contributions to A-point functions in the squeezed limit are ordered by the ratio of wavenum¬ 
bers ki/k, which is assumed to be much less than one. The lowest order contributions, up 
oc (fej/fe) _1 when the A-point function is written in terms of the overdensity <5, are fixed by 
the requirement that a uniform potential perturbation as well as a uniform velocity (boost) do 
not lead to any locally observable effect on the density held, as demanded by the equivalence 
principle [3, 4, 9, 16]. They are also referred to as “kinematical contributions”. Here we focus 
on the next order contribution, oc ( ki/k )°, which is the lowest order at which a physical cou¬ 
pling of long- and short-wavelength modes happens. More precisely, the contributions at this 
order correspond to the impact of a uniform long-wavelength density or tidal perturbation. 
When considering equal-time A-point functions, which we do throughout, and subhorizon 
perturbations ki 3> aH , the kinematical contributions disappear, and the physical ( ki/k)° 
contributions are the leading contribution to the A-point function in the squeezed limit. 1 

1 This extends to ki < aH if the density perturbation is written in synchronous-comoving gauge [16]. 
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In this paper, we disregard tidal fields, which leads us to first angle-average over the 
IV —2 small momenta (wavenumbers) in the IV-point function. Specifically, we consider Sn -2 
defined through 

S N -2(k, k'- h, • • • , k N . 2 ) = Cl ~ k ^ 1 <<y(k)(5(k , )<5(k 1 ) • • • 5(kw_ 2 )>' , (1-1) 


where k, are unit vectors and (d(ki) • • • <5(kjv))[. denotes the nonlinear connected matter IV- 
point function with the momentum constraint (27r) 3 d'£)(ki + ■ ■ ■ + kjv) dropped. Note that 
the momentum constraint fixes k' in terms of k and ki,..., k n- 2 - We now let k\,..., kjst- 2 
go to zero, and normalize the result by the nonlinear power spectrum P(k ) and the linear 
power spectra Pi(k \) • • • Pi(kN- 2 ) to obtain a dimensionless quantity: 


RN-2(k) 


lim 

hi —>-0 


S N - 2 {k,k']kir ■ -k N -2) 
P(k)P l (k l )---P l {k N .2) ' 


( 1 . 2 ) 


Note that in this limit, spatial homogeneity enforces k' = —k + 0(kj./k), so that (for statisti¬ 
cally isotropic initial conditions) the r.h.s. only depends on k. In Appendix A (see also [10]), 
we show that the R n (k) exactly correspond to the power spectrum response functions, which 
quantify the change in the nonlinear matter power spectrum to an infinite-wavelength density 
perturbation. These response functions are defined as the coefficients of the expansion of the 
power spectrum in the linearly extrapolated initial overdensity 5 lo'- 


P(k,t\5 L0 ) 


E l r * in 

—R n (k, t) 5 L0 D(t) P(k,t), 

n\ L J 

n=0 


(1.3) 


where P{k,t\5u)) is the nonlinear matter power spectrum at time t in the presence of a 
homogeneous (infinite-wavelength) density perturbation, and D{t) is the linear growth factor 
normalized to unity today. We have set Ro(k,t) = 1 by definition. Thus, by measuring R n , 
we measure the angle-averaged squeezed limit (Eq. (1.2)) of the nonlinear matter (n+2)-point 
function. 

For n = 1, the response R\ describes the angle-averaged squeezed limit bispectrum. 
This relation has been derived several times in the literature (e.g., [11, 14, 17]). Refs. [18, 19] 
considered the case of n = 1 in the context of the power spectrum covariance. Ref. [10] 
considered the general angle-averaged case of n long modes and l short modes; however, a 
hierarchy between all of the long modes ki <C ki+± was assumed which we do not assume 
here. 

Independently of the derivation of Eqs. (1.2) and (1.3), we here present accurate mea¬ 
surements of R n for n = 1, 2, 3 using N-body simulations which do not rely on approximations. 
Specifically, we resort to N-body simulations with an external homogeneous overdensity im¬ 
posed via the separate universe approach described in Ref. [20] (see also [19, 21, 22] when 
the overdensity can be approximated to be small). A flat FLRW universe with a homo¬ 
geneous overdensity is exactly equivalent to a different, curved FLRW universe [16, 23], so 
that N-body simulations in this modified cosmology provide, in principle, the exact result 
for the response functions R n (k). This in turn corresponds to the exact (in the limit of infi¬ 
nite volume and resolution) measurement of the squeezed-limit IV-point function (Eq. (1.2)). 
Ref. [19] presented simulation measurements of R\{k). In Ref. [20], in which we introduced 
the separate universe simulation technique, we already briefly presented a subset of the re¬ 
sults shown here (the so-called “growth-only response” defined later in this paper) for n = 2 
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and 3 as a sample application of the method. Here, we present the full results for n = 2 and 
3 for the first time. 

Many semi-analytical approaches to nonlinear large-scale structure assume that non¬ 
linear matter statistics can be described as a unique function of the linear matter power 
spectrum, i.e. the power spectrum of initial fluctuations linearly extrapolated to a given 
time. In the context of consistency relations, this approximation has been studied in, e.g., 
Refs. [10, 11]. This ansatz is motivated by the fact that in Einstein-de Sitter (flat matter- 
dominated universe), and to a very good approximation in ACDM, the perturbation theory 
predictions factorize into powers of the linear growth factor and convolutions of products 
of the initial matter power spectra and time-independent functions. Another way to phrase 
this ansatz is that nonlinear large-scale structure only depends on the normalization of the 
fluctuations at a given time, and not on the growth history. In the context of squeezed-limit 
Appoint functions, this ansatz can be tested quantitatively by comparing the outputs of sep¬ 
arate universe simulations at a given time with simulations in which the initial amplitude of 
fluctuations is rescaled to match the linear power spectrum at the same time. The difference 
between these “rescaled initial amplitude” simulations and the separate universe simulations 
corresponds to the error made in the ansatz of assuming that the linear power spectrum at 
a given time uniquely describes nonlinear large-scale structure at the same time. Ref. [24] 
studied this for n = 1 and found that the two simulations differ in the nonlinear regime. 
Ref. [13] performed a closely related test using the matter bispectrum. In this paper, we 
study this comparison in more detail and for n = 1,2 and 3. 

The outline of the paper is as follows. We develop semi-analytic predictions for the 
power spectrum response in section 2. We then describe the N-body simulations used in this 
paper in section 3. Results and comparisons are presented in section 4, and we conclude in 
section 5. The appendices present the proof of eqs. (1.2) and (1.3) (Appendix A) as well as 
various useful results on the separate universe picture necessary for the analytical approaches 
in section 2. 


2 Power spectrum response 


We define the nth-order response function R n {k ) of the power spectrum as the nth derivative 
of the power spectrum with respect to the linearly extrapolated (or Lagrangian) overdensity 
5l, normalized by the power spectrum. The definition consistent with Eq. (1.3) is 


Rn(k,t) 


1 d n P(k,t\S L ) 
P(k ) d(6 L (t))" 


Sl=0 


( 2 . 1 ) 


where 5i,(t) = 5LoD(t). In the following, we will frequently suppress the time argument for 
clarity. Analogously, one can define the power spectrum response functions with respect to the 
fully evolved (or Eulerian) nonlinear overdensity 5 p . Since we can expand the nonlinear over¬ 
density in powers of 5 l with known coefficients via the spherical collapse (see Appendix B.l), 
the nth-order Eulerian response function is given by a sum of R m with m < n. In this paper, 
motivated by the relation Eq. (1.2), we mainly consider the Lagrangian response functions. 
In the remainder of this section, we develop semi-analytic models for the response functions 
based on the separate universe picture. 
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2.1 Separate universe picture 

An infinite-wavelength adiabatic density perturbation 5 P behaves like an independent curved 
“separate” universe [21-23, 25, 26], in which 5 P is absorbed in a modified background matter 
density by a modification of the cosmological parameters. A positive overdensity implies a 
slower expansion, which can be quantified by the relative difference in the scale factor of 
the background and modified cosmology 5 a (t) = a(t)/a(t) — 1, where here and throughout 
the paper a tilde denotes quantities in the modified, separate universe cosmology. On the 
other hand, a(t) refers to the fiducial cosmology, for which we aim to calculate the response. 
Consequently, the comoving coordinates of the two cosmologies are related by 

x= a ^5t=[l + d a (t)]5t. (2.2) 

a{t) 

Furthermore, due to mass conservation, the fractional difference in the scale factor is related 
to the overdensity 6 P by 

1 + 5 p (t) = [1 + 5 a (t)] 3 . (2-3) 

Using the separate universe picture, we can regard the matter power spectrum in this patch 
just as that of a region with no homogeneous overdensity but properly modified cosmology. 
The modification of the cosmology is such that the shape of the linear power spectrum is 
unchanged, since the ratio of photon, baryon, and cold dark matter densities is unmodified; 
moreover, the transfer function parameters are unchanged: Cl m h 2 = Q m h 2 and £li,h 2 = H;,/? 2 . 
Thus, only the growth of structure is affected. 

The power spectrum that enters in the response given by Eq. (2.1) is defined with respect 
to the background density and comoving coordinates of the fiducial cosmology. Hence, the 
power spectrum calculated for the modified cosmology has to be mapped to that with respect 
to the background density and comoving coordinates of the fiducial cosmology. This mapping, 
described in more detail below, yields the “reference density” and “dilation” contributions to 
the response [11, 14, 17, 24], These can be calculated exactly at any scale k to any given order 
given the nonlinear matter power spectrum in the fiducial cosmology. That is, we do not 
need to run separate simulations to calculate these effects. They are thus merely “projection 
effects”, unlike the effect of the modified cosmology on the growth of structure, which requires 
a simulation in order to provide an accurate estimate. Let us denote the power spectrum for 
the modified cosmology as P(k). Then, the reference density effect simply rescales the power 
spectrum, 

p(fc) ref. density [1 + ^ ]2j p (fc)) (2 . 4) 

where the argument of P(k) is not modified. The dilation effect due to the change in the 
coordinates given by Eq. (2.2) implies k —> k = (1 + 5 a )k and changes the power spectrum 
by (see Appendix D) 

P(k) dil = ion [1 + <5 a ] 3 P([ 1 + Sa]k ). (2.5) 

Putting the two together and using Eq. (2.3) yields 

P(k) = [l + S p }P{[l + 6 a \k) , (2.6) 

where all quantities are evaluated at some fixed time t. Note that one prefactor of 1 + 5 P 
cancels, since the effect of the increased density is partially canceled by the corresponding 
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decrease in physical volume. For a flat matter-dominated fiducial cosmology, it is straight¬ 
forward to derive series solutions for 5 P and 5 a (Appendix B.l) of the form 


OO 

^a(^) = ^ ^ 

n =1 


^LoD(t) 


n 


Sp(t)=J2fn 

n =1 


<■$LoD(t ) 



(2.7) 


where D(t ) = D(t)/D(to) is the fiducial growth factor normalized to one at the epoch to to 
which we extrapolate 5lo = £l(^o)> and e n , f n are rational numbers. 

The third contribution to R n comes from the effect of the modified cosmology on the 
growth of structure, which as mentioned above is the physical contribution which requires N- 
body simulations for an accurate measurement. We thus define a set of growth-only response 
functions G n (k) which isolate the nontrivial effect of the long-wavelength perturbation on 
the growth of small-scale structure, 


G n (k) 


1 d n P{k) 
P(k) dd n L 


6 l =o 


( 2 . 8 ) 


That is, G n are defined as R n without the contributions from the reference density and 
dilation given by Eq. (2.6). This definition is an extension of the similar decomposition for 
n = 1 shown in refs. [11, 14, 17, 24], Thus, the formula for the power spectrum (w.r.t. global 
coordinates) in the presence of a long-wavelength overdensity is given by 


P(k\S L ) 


[1 + d p ] 


1 + E fp G n(W n L P(k) 


n= 1 


k=[l+6 a ]k 


(2.9) 


Clearly, by the Leibniz rule, at any given order n the total or “full” response R n (k ) [Eq. (2.1)] 
is composed of the functions G m (k) and the numbers e m , f m with 1 < m < n, where the 
e m multiply derivatives of Gi(k) and P(k) with respect to k (up to the n-th derivative). 
Specifically, the first three full response functions are given by 


Ri(k) 

R2(k) 

2 


Rs(k) 

6 


kP'(k) . . 

fl + ei ^w + l( 

f 2 + + e? + GM + hei kP\k) 


( 2 . 10 ) 


+ fiGi(k) + e\ 


P(k) ' 2 P(k) 

kP'(k) 


P(k ) 


m 


Gi(k) + e 1 kG' 1 (k), 


( 2 . 11 ) 


kP'(k) , , , G 3 (k) , kP\k) , , G 2 (k) , , kP'{k) 

hGl[k)ei TW + h + “6“ +63 W ++/l62 ~m 

+ +hGl(k) + h p_m + (/iei+ e2)kG>l[k)+ 


+ 


2 P(k) J J P(k) 

G' 2 (k) 2 kP'(k),^, ... 3 k 3 P"'(k) 0 k 2 P"(k) 

-WW + 2eie2 w 


+ ei 


kP'(k) G 2 [k ) 
P{k) 2 


+ G 1 (k)(e 2 


kP'(k ) 2 k 2 P”{k) 
P(k) +£l 2 P(k) 


( 2 . 12 ) 


where the primes denote derivatives with respect to k. 
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2.2 Linear power spectrum predictions 

We now evaluate Eq. (2.9) for the simplest case, i.e., the response of the linear matter power 
spectrum. In linear theory, the growth is scale-independent and given by the linear growth 
factor. Thus, the growth-only response functions are scale-independent and just described 
by the linear growth factor in the modified cosmology D(t), 


/'-(linear = 1 <**(&) 

D 2 d5l 


(2.13) 


A perturbative expansion of D in powers of 5l for a flat matter-dominated fiducial cosmology 
is derived in Appendix C, with the result given in Eq. (C.18), 


D(t) = D(t) 


OO 

1 + ^2 9 n 

n= 1 


Slo D(t) 


(2.14) 


Thus, for an Einstein-de Sitter fiducial universe (and to high accuracy in ACDM), the linear 
response functions are simply constants. Inserting the result from Eq. (C.17), we obtain 


{ 


G 


linear 

n 


} 


n=l,*”4 


26 3002 240272 197919160) 
21’ 1323’ 43659 ’ 11918907 J ' 


(2.15) 


Eq. (2.9) evaluated for the linear matter power spectrum Pi(k,t ) then becomes 


Pi(k,t\S L ) 


[1 + 5 p (t)] 



2 

-Pz,fid([l + S a (t)]k,t). 


(2.16) 


Inserting the series expansions derived in Appendix B and Appendix C, we obtain 

/ OO \ / OO „\ ^ 

Pi(k,t\S L0 ) = I 1 + ^2 f n [6 L oD{t)] n J (l + ^2g n d L0 D(t) 


\ 71=1 

x -Pi,fid ( 


71=1 


1 + J2 e ^L 0 D(t)Y 


71=1 


k, t 


(2.17) 


Eq. (2.17) allows for a consistent expansion in 5lq. Specifically, d n Pi(k )/is given by the 
n-th order coefficient in this expansion, multiplied by n!. 


2.3 Nonlinear power spectrum predictions 

Beyond the linear matter power spectrum, the growth coefficients G n will become scale- 
dependent functions G n (k). Consider now what standard perturbation theory (SPT) pre¬ 
dicts. The power spectrum prediction is given by a series 

P spt (A:) = P t (k) + P 1 - loop (fe) + P 2 - loop {k) + • • • , (2.18) 

where P n_loop scales as [Pi] n ■ In an Einstein-de Sitter universe, one can show (e.g., [27]) that 
the time- and scale-dependence of each order in perturbation theory factorizes, so that one 
can write 

P SPT {k , t) = b 2 {t)Pt{k , t 0 ) + b^tjP 1 - 100 *’^, t 0 ) + b 6 P 2 - lo °P{k, t 0 ) H-, (2.19) 


- 7 - 














where P n ~ loop (k, to) is a convolution of n factors of Pi{k,to) with time-independent coeffi¬ 
cients. While Eq. (2.19) is only strictly correct in Einstein-de Sitter, it is used very commonly 
for ACDM as well, since departures from the exact result are typically of order 1% or less, 
and since it simplifies the calculation significantly. Various variants of SPT, such as the 
renormalized perturbation theory (RPT) [28], share the same property. 

In the context of this paper, Eq. (2.19) allows for a very simple evaluation of the growth- 
only response: as discussed above, the shape of the linear power spectrum in the modified 
cosmology is unchanged, and hence P SPT (fe) can be simply evaluated by replacing the fiducial 
D(t ) in Eq. (2.19) with the modified one, Eq. (2.14). This is equivalent to assuming that 
the entire late-time cosmology dependence of the nonlinear matter power spectrum enters 
through the linear growth factor [10, 11, 14]. 

Apart from the SPT calculation, we can also apply this approximation to any prescrip¬ 
tion that maps a given linear power spectrum to a nonlinear one. In particular, we will show 
results for halofit [29]. In this case, where the dependence on the linear growth factor is 
not explicit, we instead compute the derivative with respect to the normalization of the linear 
power spectrum, 


d das d 
dD dD das ’ 


( 2 . 20 ) 


which at the redshift considered yields the equivalent change of the linear matter power 
spectrum. This leads to 


d n P{k) 

dD n 


a 8 


, d n P(k ) 

da £ 


( 2 . 21 ) 


We use a five-point stencil with a step size of 0.75% in as to compute numerically the 
derivatives with respect to as- In conjunction with the change of the linear growth factor 
Eq. (2.14), this allows us to compute the growth-only response G n (k ) for perturbation theory 
as well as fitting formulae of the nonlinear matter power spectrum. 

Further, we can test this prescription to all orders in SPT calculations, and indepen¬ 
dently of fitting functions, by performing simulations with a rescaled initial power spectrum. 
This is the subject of section 3.2. 


2.4 Halo model predictions 

In the halo model (see [30] for a review), all matter is assumed to be contained within 
halos with a certain distribution of mass given by the mass function, and a certain density 
profile. Along with the clustering properties of the halos, these quantities then determine 
the statistics of the matter density field on all scales including the nonlinear regime. IV-point 
functions can be conveniently decomposed into 1- through V-halo pieces. In the following, 
we will follow the most common halo model approach and assume a linear local bias of the 
halos. This is the most popular choice in the literature, although it can clearly be improved 
upon (and is not strictly consistent, as we will see). 

Adopting the notation of Ref. [18], the halo model power spectrum, PuM(k), is given by 

Pu M (k) = P 2h (k) + P lh (k), (2.22) 

P 2h (k)= [Il(k)] 2 m ), 

P lh (k)=I°(k,k), 
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where 


!%{k i, • • • k m ) = / d In M n(lnM) 


M 

P 


(2.23) 


b n (M ) u(M\k\) ■ ■ ■ u{M\k n 


and n(lnM) is the mass function (comoving number density per interval in log mass), M is 
the halo mass, b n {M) is the n-th order local bias parameter, and u(M\k ) is the dimensionless 
Fourier transform of the halo density profile, for which we use the NFW profile [31]. We 
normalize u so that u(M\k —> 0) = 1. The notation given in Eq. (2.23) assumes &o = 1- 
u(M\k) depends on M through the scale radius r s , which in turn is given through the mass- 
concentration relation. All functions of M in Eq. (2.23) are also functions of 2 although we 
have not shown this for clarity. In the following, we adopt the Sheth-Tormen mass function 
[32] with the corresponding peak-background split bias, and the mass-concentration relation 
of Ref. [33] . The exact choice of the latter only has a small impact on the predictions which 
does not affect our conclusions. A dependence of halo profiles on the background cosmology 
can however change the response functions on small scales k > 1 /i Mpc^ 1 (see section 2.4.2 
below). The derivations of the halo model response given below generalize the linear response 
calculations of [18, 34] to arbitrary nonlinear order. 


2.4.1 Total halo model response 

We now derive how the power spectrum given in Eq. (2.22) responds to a homogeneous 
(infinitely long-wavelength) density perturbation 5l- For this, we consider the 1-halo and 2- 
halo terms separately. The key physical assumption we make is that halo profiles in physical 
coordinates are unchanged by 5l- That is, halos at a given mass M in the presence of 5l 
have the same scale radius r s and scale density p(r s ) as in the fiducial cosmology. We will 
discuss this assumption in section 2.4.2. Given this assumption, the density perturbation 
5l then mainly affects the linear power spectrum, which determines the halo-halo clustering 
(2-halo term), and the abundance of halos at a given mass. 

We begin with the 2-halo term. The response of the linear power spectrum was derived 
in Eq. (2.17) in the previous section. The expression for the 2-halo term in Eq. (2.22) is simply 
the convolution (in real space) of the halo correlation function in the linear bias model with 
the halo density profiles. By assumption, the density profiles do not change, hence l\ only 
changes through the bias b\(M) and the mass function n(lnM). The bias bw(M) quantifies 
the A-tli order response of the mass function n(lnM) to 5 l [35, 36]: 


b N (M) 


1 8 N n(\nM) 
n(lnM) 85% 0 ’ 


so that 


8 N n(\nM ) 
85% o 


bisr(M)n(\n M). 


Thus, 


d N 

ds" 


I\(k) 


s L =o 


j dlnM 



p)N 

-^MM)nQnM)] 


u(M\k ) 

S L =0 


I% + \k). 


(2.24) 


(2.25) 


Note that in the large-scale limit, k —> 0, this vanishes for N > 1 by way of the halo model 
consistency relation 


J d In M n(ln M) 0Q b N (M) = j J’ * > j ; 
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(2.27) 


For finite k however, Eq. (2.25) does not vanish. We thus have 

OO . 

ll(k,t\S L0 ) = J2 -J[ l+1 (k,t)[D(t)6 L0 ] n . 

ZJ n! 

n =0 

Thus, the two-halo term in the presence of 5lq becomes 


P 2h (k,t\8 L0 ) 


1 + E fn[^L 0 D(t)] n 1 + £ g n \5 L0 D(t ) 


71=1 
OO 


71=1 
2 


X E^ + HM)[£(0M n Pi,fid 


\n =0 


l + J2 e ^L 0 D(t)Y 


71=1 


(2.28) 



Note that we recover the tree-level result given in Eq. (2.17) in the large-scale limit. Strictly 
speaking, this expression is not consistent, since the term if implies a non-zero 62 while in 
Eq. (2.22) we have assumed a pure linear bias. Note that in Eq. (2.28) the dilation effect only 
enters in the linear , not 2-halo, power spectrum. This is a consequence of our assumption 
that halo profiles do not change due to the long-wavelength density perturbation. 

We now turn to the one-halo term. Given our assumption about density profiles, this 
term is much simpler. The only effect is the change in the mass function, which through 
Eq. (2.24) becomes 

d N 

—I»( k ,k) = I?(k,k). (2.29) 

We thus obtain 

OO . 

P lh (k,t\d L0 ) = £ -l2(k,k,t)[D(t)6 L 0 ] n . (2.30) 

71=0 


Putting everything together, we obtain 

/ OO \ / OO 

P nM (k,t\6 L 0 ) = 1 + y/„MW] n 1 + X>n \SwD(t) 


71=1 

OO 


71=1 
\ 2 


x £ S t“ + 1 (M)[»(().5Lo]” Pi, 


fid 


Vn=0 
00 


1 + e n [&Lf)D(t)Y 


n=L 


k, t 


+ '£- I 2(k,k,t)[D(t)6 LO ]' 
n\ 


(2.31) 


71=0 


The contribution oc If +1 (for n > 0) is numerically much smaller than the other terms (see 
also the discussion in section 4.2.4 of [34]). Since it is much smaller than the overall accuracy 
of the halo model description, we will neglect it in the following. This yields 


P UM (k,t\5 L0 ) 


1 + J2fn[S L oD(t)] n ) 1 + X> n \SuoD(t) 


71=1 


71=1 


x (ll(k,t)) 2 P lM 


OO 

1 + e n [5LoD(t)]' 1 

71=1 


+ ^2^ I 2(k,k,t)[D(t)S L0 ] n . 

n =0 



(2.32) 
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Explicitly, the first and second order full response functions are given by 


Rf M {k ) 

Rf M (k) 


fi + 2gi + ei 


dlnPi(k, t ) 
din k 


P 2h (k,t) + I%(k,k,t) 


2/2 + 2 /ipi + (/1 + 2 pi)ei 


dlnP ; (/c, t) 
din k 


d]nPi(k,t) 2 1 d 2 Pi(k,t ) 

2 d In k l P d(lnfc) 2 


->2h 


+ 2 .<7 2 + 4^2 
(k,t) + I 2 {k,k,t). 


(2.33) 


2.4.2 Growth-only response 

We also derive the growth-only response functions in the halo model approach. Since the 
halo profiles are assumed fixed in physical coordinates, this means that we need to rescale 
the halo model terms, 1^, accordingly. Following our discussion in section 2.1, we have 
k = (1 + 5 a )k , where k is the comoving wavenumber with respect to the modified cosmology. 
We then obtain 


growth only 


(&!>■■* ^m) — I-n 


h 


k r 


physical \ 1 + 6 a (t) ’ 1 + 6 a (t) 


(2.34) 


Inserting this into Eq. (2.31) and performing a series expansion of 5 a in Sl then allows us to 
derive the growth-only response functions G^ M (k). Note that the NFW profile we assume 
is uniquely determined by the scale radius r s (M ) for a halo of mass M, which enters the 
coefficients defined in Eq. (2.23) in the combination kr s (M). Thus, it is easily possible to 
include a dependence of the scale radius r s (M), or equivalently the halo concentration, on 
the long-wavelength density in a similar way. We will leave this for future work. 

Quantitatively, the main contribution of the rescaling Eq. (2.34) is from the 1-halo 
term oc fc), i.e. the term in the last line of Eq. (2.32). The rescaling of the other 

instances of 1^ only changes the response at the sub-percent level and we will neglect them 
in the following. We then obtain for the growth-only contribution to the halo model power 
spectrum 


iHM/ 


(k, t\5 L0 ) gr ° w = ° nly 1 + ^ g n [<5 L0 ^(t)]" {A 1 [A, t] }“ Pi M (k, t) 


n= 1 


- 

+ E *) ^o, t) k, t] [D(t)5 L0 ] n , 

n =0 


(2.35) 


where 

( OO 

1 + Y J e n[dwD(t)] n 

n— 1 

This completes the derivation of the halo model response functions. 



(2.36) 


3 N-body simulations 

Before describing the separate universe simulations in detail, we summarize features common 
to all. All simulations are gravity-only simulations and are carried out with Gadget-2 [37]. 
The starting redshift is z = 49 and the initial displacement field is computed using second- 
order Lagrangian perturbation theory. For each simulation, the particle load is 512 3 . For the 
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fiducial cosmology (5lo = 0), we choose a flat ACDM cosmology with cosmological parameters 
consistent with the current observational constraints: Q m = 0.27, h = 0.7, Qf,h 2 = 0.023, 
n s = 0.95, us = 0.8, and a comoving box size of 500 /i -1 Mpc. 


3.1 Separate universe simulations 

Using the separate universe approach presented in Ref. [20], we simulate separate universes 
corresponding to the linearly-evolved present-day overdensities of 5lo = 0, ±0.01, ±0.02, 
±0.05, ±0.07, ±0.1, ±0.2, ±0.5, ±0.7, and ±1. Then, for the separate universes, the Hubble 
constant and the curvature fraction vary between h: 0.447 to 0.883 and Qk- —2.45 to 0.372, 
respectively. The physical densities (l m h 2 , Cl\h 2 , and Ct^h 2 as well as n s and the amplitude 
of the primordial curvature power spectrum remain the same. 

The initial conditions are set up as described in Ref. [20]. For each overdensity Slo, 
we run the same 16 realizations of the Gaussian initial density field. Hence, by comparing 
relative differences between different 5lq values but the same realization, most of the sample 
variance cancels out. The 16 realizations allow us to estimate the residual statistical error. 

Given a fixed box size for the fiducial cosmology, there are two reasonable choices for the 
box sizes of the modified cosmologies. Either we match the respective comoving box sizes, i.e. 
the box size is 500 h/h in units of h _1 Mpc comoving, or we choose the box sizes such that their 
physical sizes coincide with that of the fiducial simulation at one specific output time t Qut , 
i.e. 500 ha(t out )/[ha(t ou t)} in units of h -1 Mpc comoving, where a and a are the scale factors 
of the fiducial and modified cosmology, respectively. The former choice is adequate if we are 
interested in the power spectrum response functions at the same comoving wavenumber, i.e. 
without the “dilation” effect. By using the mean density of the separate universe cosmology 
as the reference density when computing the power spectrum, we are further removing the 
“reference density” effect and are left with the growth-only response. In Ref. [20], we have 
run simulations with this choice of box-size-matching to measure the growth-only response 
functions. Here, we reproduce the simulations with a higher mass resolution and compare 
the results with the models presented in this paper in section 4.1. 

In order to measure the full response functions, we run simulations for which we match 
the physical box size. We focus on two different output times t ou t corresponding to z = 0 
and 2 : = 2 in the fiducial cosmology. As the physical size can only be matched at one specific 
time, we have to run a new set of simulations for each output time. The results of these 
simulations are presented in section 4.3. 


3.2 Simulations with rescaled initial amplitude 

We also investigate how well the effect of a homogeneous overdensity on the growth of struc¬ 
ture can be modelled by a change in the amplitude of the linear power spectrum. To this end, 
we additionally run a set of simulations for which we always assume the fiducial cosmology 
but vary the amplitude of the initial power spectrum. Specifically, for each 5lo value for 
which we simulate a separate universe, we also simulate the fiducial cosmology with the ini¬ 
tial power spectrum amplitude multiplied by D(to) 2 / D(to ) 2 , where D(to ) is the linear growth 
factor in the corresponding separate universe cosmology. The results of these simulations are 
shown in section 4.2. 
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4 Results 


For the power spectrum computation, we first estimate the density contrast S(x) on a 1024 3 
grid using the cloud-in-cell mass assignment scheme, then apply a Fast Fourier transform, 
and angular average the squared amplitude |<5k| 2 - The density contrast 5(x) = p{x)/p — 1 
describes the overdensity with respect to the reference density p. When we are interested 
in the growth-only response function, p is equal to the mean density of the separate uni¬ 
verse. When we compute the full response function, p is equal to the mean density of the 
fiducial cosmology. Similarly, for the growth-only response, distances are measured using the 
comoving coordinates of the respective cosmology 2 , whereas, for the full response, the power 
spectrum is always measured in comoving coordinates of the fiducial cosmology. 

We only report results up to a maximum wavenumber of 2/i -1 Mpc. A convergence 
study with simulations with 8 times lower mass resolution shows differences in G\ . G2 and 
G 3 of only 1 (3) to 5 (10) percent at z = 0 (z = 2) up to that wavenumber, where the 
deviations increase from the linear response function to the higher-order response functions. 
The results for the full response functions R\, R2 and R3 are converged to an even better 
degree. We therefore expect that the simulation results presented in this paper are converged 
to a sub-percent to a few percent level. 

In order to compute the first three response functions, we fit a polynomial in Sl to the 
fractional difference in the measured power spectrum A k{^>L) = P(k\5i,)/P{k\$L = 0) — 1 
for each £;-bin. For the fit, we only include results from separate universe simulations with 
|<5i(tout)| < 0.5 and use a polynomial with degree 6 to be unbiased from higher-order response 
functions. As the random realization of the initial density field is the same across different 
5 l values, the corresponding power spectra are strongly correlated. By considering the ratio, 
or the relative difference, of two power spectra a large fraction of the noise cancels. However, 
for the same realization, the measured fractional differences A^(^l) are still correlated over 
5 l ■ As the number of realizations (16) is not large enough to reliably estimate the covariance 
between different Sl values, we cannot include this correlation in the polynomial fitting. In¬ 
stead, we construct quasi-decorrelated samples of Afc(<5 l) by randomly choosing a realization 
for each Sl value. Fitting many of those subsamples allows for a robust error estimation of 
the derived response functions. 

4.1 Growth-only response functions 

Figure 2 shows the first three growth-only response functions measured from the simulations 
at z = 0 (left column) and z = 2 (right column). These correspond to the fully nonlin¬ 
ear squeezed limit bispectrum (3-point function), trispectrum (4-point function) and 5-point 
function, and are essentially the same as in figure 2 of Ref. [20], although the results pre¬ 
sented there are derived from lower resolution simulations (particle load of 256 3 ). The small 
wiggles in the growth-only response functions result from the damping of the baryon acoustic 
oscillations (BAO), which depends on the amplitude of density fluctuations and thus on Sl- 

Let us compare the simulation results to the theoretical predictions discussed in sec¬ 
tion 2. On sufficiently large scales, the perturbation theory predictions are the most accurate, 
as expected. At high redshift, the 1-loop predictions best describe the results overall. The 
1-loop predictions also show a BAO damping effect. At z = 0, the growth-only response 
is captured best by the halofit prescription (in case of Gi) or the halo model (in case of 
G2, G3). We see that the halofit prescription describes the simulation results of the linear 

“Note, however, that the unit of length is always /i _1 Mpc, where h corresponds to the fiducial cosmology. 
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Figure 2: The first three growth-only response functions of the power spectrum measured 
from the separate universe simulations at z = 0 (left) and z = 2 (right). The error bars show 
the statistical error derived by random resampling of the data (see text). For data points 
apparently without error bars, the statistical error is smaller than the size of a dot. 


response well at both redshifts, but performs significantly worse for the higher-order response 
functions. The BAO damping effect is essentially absent in both halofit and halo model 
predictions. Overall, none of the models is able to accurately describe the simulation data 
in the nonlinear regime, with discrepancies at z = 0 ranging from 20% in the best case to a 
factor of several. These discrepancies are not surprising given that we are looking at scales 
beyond the validity of perturbation theory and at higher Appoint functions for which the 
senri-analytical approaches were not tuned. 

The halo model prediction does not asymptote exactly to the linear result in the k —> 0 
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limit. This is because the 1-halo term asymptotes to a white noise contribution in this limit, 
and since the 1 -halo term contributes to G n due to the dependence of the halo mass function 
on 5l (section 2.4), this induces a correction to the linear prediction which contributes on 
large scales. Physically, this occurs because the halo model does not enforce momentum 
conservation of the matter density field. This issue can be fixed by introducing a “mass 
compensation scale” [38]. 

The halo model predictions can be tuned to better match the simulation results by al¬ 
lowing for a dependence of the halo profiles on the long-wavelength density, which is expected 
on physical grounds (see also [24] ). Specifically, if the the scale radius of halos at fixed mass 
increases in the presence of a long-wavelength density perturbation, this lowers the peak in 
the response and thus could lead to better agreement with the simulations results. A detailed 
investigation of this is beyond the scope of the present paper. 

4.2 Comparison to simulations with rescaled initial amplitude 

All models for the growth-only response functions that we have presented in section 2 and 
shown above are based on the approximation that we can trade the effect of 5l for an 
appropriate change to the linear growth factor (or equivalently, the linear power spectrum). 
But how well does this approximation work? Using the set of simulations described in 
section 3.2, we can explicitly test this approximation on all scales including the nonlinear 
regime. 

In figure 3, we show the growth-only response functions measured from two different 
sets of simulations. In case of G i, this comparison was also shown in figure 6 of [24], and our 
results agree with theirs . 3 The “rescaled amplitude simulations” of section 3.2 all assume 
the fiducial cosmology but vary the amplitude of the linear power spectrum used to initialize 
the simulations so as to match the linear power spectrun in the modified cosmology at the 
given output times [Eq. (2.14) and Eq. (2.20)]. On linear scales, these simulations thus agree 
with the “separate universe” simulations by construction. As the simulations share the same 
random realization of the initial density, the sample variance (noise in the upper panels) gets 
vastly reduced when considering the difference of the measured response functions, A G n = 
G<rescaied — Gn parate . This is shown in the lower subpanels of figure 3, where we have divided 
A G n by the corresponding linear growth-only response, i.e. the prediction in the k —> 0 limit. 

The differences seen in figure 3 are caused by the different growth history, which is not 
captured by the rescaling of the initial amplitude. Following the discussion in section 2.3, the 
commonly used SPT approach factorizing the growth factor and scale dependence assumes at 
all orders that a long wavelength density perturbation enters exclusively through the modified 
linear growth. Thus, even when calculated to all orders, the best that this SPT calculation 
could do is to reproduce the rescaled amplitude result in figure 3, which deviates from the 
actual response at z = 0 by 10% at k ~ 0.5/tMpc ^ 1 and 20% at k ~ 1/iMpc ^ 1 for Gi, 
and significantly worse for the higher-order response functions. At 2 = 2 on the other hand, 
the rescaled-amplitude G\ matches the separate universe response to better than 10% even 
beyond k = 1 /iMpc -1 , and for G 2 , G 3 performs significantly better as well. 

There are two possible explanations for these discrepancies in the SPT context. First, 
using the SPT kernels derived for an Einstein-de Sitter universe (which have time-independent 
coefficients), with the ACDM linear growth factor replacing the Einstein-de Sitter a(t), could 
become highly inaccurate for ACDM at higher orders. Note that the same issue exists for a 

3 Note that the authors of ref. [13] perform a different comparison using the time derivative of the nonlinear 
power spectrum in simulations of the fiducial cosmology. 
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Figure 3: Comparison of the growth-only response functions G i, G 2 , G 3 (top to bottom) 
measured at z = 0 (left column) and z = 2 (right column) from one realization of the 
separate universe simulations and from the same realization simulated by varying the initial 
amplitude. The bottom sub-panels show the difference, A G n = G^escaied _ ^separate^ divided 
by the response of the linear matter power spectrum, Gj) near . 
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fiducial flat Einstein-de Sitter universe, since for Sl / 0 the quantity fl m // 2 is no longer 1 (in 
fact, d(Q m /f 2 )/d§L = —5/21 [14]; see also the discussion in [13]). There is no indication of 
such a strong effect at low orders in perturbation theory, where this approximation typically 
performs to better than a percent [27]. Furthermore, ref. [14] found that a cancelation in the 
curvature contribution to the growth integral suppresses this effect. Finally, ref. [24] shows 
that the growth-only response of the power spectrum to a change in the Hubble constant 
while keeping Q m h 2 fixed follows the separate universe response very closely (Fig. 6 there). 
If the much larger discrepancies between separate universe response and rescaled amplitude 
response were due to the cosmology dependence of the SPT kernels, one would not expect 
this to be the case. Nevertheless, we do not claim to be able to rigorously exclude this 
possibility. 

The other possibility, more likely in our opinion, is that the discrepancy between rescaled 
amplitude and full separate universe simulations is due to effective non-perfect fluid terms, 
such as pressure and anisotropic stress, in the dark matter fluid [39]. The effective fluid prop¬ 
erties depend on highly nonlinear small scales which are not described by the Euler-Poisson 
system. Their value can depend on the growth history (as well as the power spectrum shape) 
thus leading to a discrepancy between rescaled amplitude and separate universe simulations. 
Assuming this interpretation is correct, figure 3 explicitly shows the breakdown of SPT on 
nonlinear scales as effective pressure, anisotropic stress and sound speed need to be included. 
Separate universe simulations can be used to measure the response of these effective terms 
to a long-wavelength overdensity, which is crucial when modeling (N > 2)-point functions. 
The results shown in figure 3 are analogous to what has been found for the mass function of 
halos which is a key ingredient in the halo model description of the nonlinear matter density 
field. The mass function shows departures from being a simple function of the linear matter 
power spectrum at the 5-10% level [40, 41]. 

In an Einstein-de Sitter cosmology with scale-invariant initial power spectrum Pi(k) oc 
k n , there is only one characteristic spatial scale at any given time, which corresponds to the 
scale at which the density field becomes order 1 [42], Let us denote this wavenumber as 
k^h{t). Then, the response functions have to follow a universal function of k/k^hit), he. 
G m (k,t) = G m {k/k-^i) (keeping the index of the initial power spectrum fixed). Thus, in 
this specific case, separate universe simulations and rescaled-amplitude simulations will give 
exactly the same result when compared at fixed fc/^NL- The departures shown in figure 3 
can thus be seen as a consequence of the ACDM background and the departure from scale- 
invariance of the initial power spectrum. It would be interesting to disentangle the two 
effects, e.g. by performing separate universe simulations in ACDM with scale-invariant initial 
conditions. We leave this for future work, but point out that when plotting the differences 
shown in the lower panel of figure 3 as a function of k/k^, we still find a factor of several 
difference in the z = 0 and z = 2 results. 

Note that in Ref. [38], the authors performed a similar comparison as the one shown 
for G\ in figure 3. However, instead of studying the growth-only response as considered 
here, they included the reference density effect and found a significantly larger discrepancy 
between the two measurements considered. As discussed in section 2, however, the reference 
density effect is a simple remapping, and the growth-only response considered here is the 
proper quantity to compare for the question we are interested in, namely to measure the 
effect of the growth history. 
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Figure 4: The first three full response functions of the power spectrum measured from the 
separate universe simulations at z = 0 (left) and z = 2 (right). 


4.3 Full response functions 

We now turn to the results for the full response functions, i.e. including the “dilation” and 
“reference density” effects. The results of the simulations and the model predictions are 
shown in figure 4. The oscillations in the response functions can be traced back to the 
BAOs in the power spectrum. The BAOs propagate to the response functions primarily by 
the “dilation” effect, which yields derivatives of the power spectrum with respect to k (see 
Eqs. (2.10)-(2.12)). The 1-loop perturbation theory predictions describe the simulation re¬ 
sults accurately up to k < 0.15 /i _1 Mpc and k < 0.3 /i” 1 Mpc at z = 0 and z = 2, respectively. 
As the other theoretical models do not include the damping of the BAOs in the nonlinear 
power spectrum, they predict oscillations in the response functions which are too large. To 
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Figure 5: The first three Eulerian response functions of the power spectrum measured from 
the separate universe simulations (data points) at z = 0 (left) and z = 2 (right). The lines 
show the corresponding linear combinations of the Lagrangian response functions using the 
f n coefficients derived for the Einstein-de Sitter universe (see Eq. (4.1) and Eq. (B.18)). 


improve the accuracy of those models around the BAO scale, one would need to put in the 
BAO damping by hand. In the nonlinear regime, none of the models is able to reproduce the 
simulation data. In principle, one could build a hybrid model for the full response by com¬ 
bining an accurate prediction of the nonlinear power spectrum of the fiducial cosmology and 
the growth-only response functions G n (k ) discussed in the previous section. In this paper, 
however, we do not pursue this approach. 


4.4 Eulerian response functions 

So far, we have always considered the response to the linearly-extrapolated initial (La¬ 
grangian) overdensity 5l■ We now consider the corresponding response to the evolved nonlin¬ 
ear (Eulerian) overdensity 5 p . Using the expansion derived for the Einstein-de Sitter universe, 
Eq. (B.18), we find 


^Eulerian^) = ? 

R™*™(k) = R 2 (k)-2f 2 R 1 (k), 

if ulerian (A0 = Rs(k) - 6 f 2 R 2 {k) + 6 (2/| - h) R\{k ). (4.1) 

In figure 5, we compare the directly measured Eulerian response functions with the appro¬ 
priate linear combinations of the measured Lagrangian response functions. The agreement 
is excellent as expected, especially at high redshift at which the ACDM universe is very well 
approximated by the Einstein-de Sitter universe. 

Interestingly, the higher-order Eulerian response functions are much smaller than in the 
Lagrangian case. That is, the response of the nonlinear matter power spectrum to a uniform 
nonlinear final-time density 6 P is close to linear. This is most likely due to the fact that 
the growth-only response functions are subdominant compared to the rescaling and reference 
density contributions, especially at higher order. In this case, Eq. (2.6) implies a close to 
linear scaling with 5 p . 
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5 Conclusions 


In this paper, we employed dedicated N-body simulations using the separate universe tech¬ 
nique presented in Ref. [20] to compute the response of the nonlinear matter power spectrum 
to a homogeneous overdensity superimposed on a flat FLRW universe. The response func¬ 
tions we computed give the squeezed limits of the 3-, 4-, and 5-point functions, in which all 
but two wavenumbers are taken to be small and are angle-averaged. By virtue of the separate 
universe technique, we reach an unprecedented accuracy of these nonlinear matter IV-point 
functions. 

The response function consists of three parts: changing the reference density with re¬ 
spect to which the power spectrum is defined; rescaling of comoving coordinates; and the 
effect on the growth of structure. The former two effects can be calculated trivially, whereas 
the third one requires separate universe simulations. We have compared the simulation re¬ 
sults with analytical and semi-analytical results, in particular standard perturbation theory 
(SPT), the empirical fitting function halofit, and the halo model, finding that SPT typ¬ 
ically yields the best results at high redshifts. The fitting function and halo model, while 
qualitatively describe the trends seen in the response functions, give a poor quantitative 
description on nonlinear scales. 

A fundamental assumption of all of the analytical and semi-analytical methods used in 
this paper, including standard perturbation theory at any order, is that nonlinear matter 
statistics at a given time are given solely by the linear power spectrum at the same time, and 
do not depend on the growth history otherwise. As was done in [24] for the response function 
for n = 1, we were able to test this assumption for n = 2 and 3 quantitatively by comparing 
the separate universe simulations with simulations with a rescaled initial power spectrum 
amplitude. We find that this assumption fails at the level of 10% at k — 0.2 — 0.5/iMpc _1 
for 5- to 3-point functions at z = 0. The failure occurs at higher wavenumbers at z = 2. In 
the context of SPT, this may signal a breakdown of the perfect fluid description of the dark 
matter density field at and beyond these wavenumbers. In other words, even if computed to 
all orders, SPT (and its variants such as RPT [28]) fails to describe the nonlinear structure 
formation beyond these wavenumbers. Therefore, our results yields a quantitative estimate 
for the scales at which effective fluid corrections become important in the bispectrum and 
higher Appoint functions, and at which one should stop trusting pure SPT calculations. 

Finally, we point out that the approach presented here can be augmented to measure 
more general squeezed-limit IV-point functions, by including the response to long-wavelength 
tidal fields and by considering the response of small-scale n-point functions in addition to 
the small-scale power spectrum considered here. 
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A Squeezed limit IV-point functions and power spectrum response 

In this appendix, we prove the relation Eq. (1.2) between the power spectrum response and 
the squeezed limit IV-point functions. We only consider equal-time IV-point functions, to 
which there are no boost-type contributions from kinematical consistency relations. Further, 
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we assume that the long-wavelength modes are well inside the horizon, removing gauge- 
dependent terms present for horizon-scale modes. Note that the relations derived here retain 
their formal validity even for horizon-scale long-wavelength modes if the matter density per¬ 
turbation is evaluated in synchronous-comoving gauge [16]. 

As in Eq. (1.3), we expand the power spectrum as a function of the linearly extrapolated 
initial overdensity 5 lo as 


P(k,t\S L0 ) 


E i , in 

—Rn(k,t) S L0 D(t) P(k,t ), 

n=0 n! L J 


(A.l) 


where R n (k,t) are response functions with Ro(k,t) = 1. At the same order in derivatives, 
that is at the same order in ki/k of the squeezed-limit Appoint function, the power spectrum 
will also depend on the long-wavelength tidal field which can be parametrized through 

i^(k)= (^-I^«5(k). (A.2) 


Exploring the response of the power spectrum to a long-wavelength tidal field is beyond the 
scope of this paper. We remove the dependence on Kij by performing an angle-average of 
the long-wavelength modes which cancels the tidal field contributions. 

In the following, we will suppress the time argument for clarity. We let S n be defined 
as in Eq. (1.1), 

S n {k,k'-k u --- ,K)= J^-J ^^(k^kO^kx)---^))' . (A.3) 


Here, the prime denotes that the factor (2-7r) 3 (5D(k + k' + ki... n ) 
Yli =l kj- We consider the limit 


lim 

hi —^0 


S n (k,k ,k\,--- fc n ) 
P{k)P l [k l )---P l {k n ) ’ 


is dropped, where ki... n = 


(A.4) 


which means that all |k*| are taken to zero. In this limit, spatial homogeneity enforces 
k' = —k + 0{ki/k ), so that (for statistically isotropic initial conditions) the r.h.s. only 
depends on k. 

In order to prove Eq. (1.2), we first note that since we are interested in the limit ki —> 0, 
we can replace <5(k ? ;) in Eq. (1.1) with the linear density field <5r,(k,;). We further transform 
kj into real space, writing 


(fc; ki, ■ ■ ■ , k n ) = J d ^r ■■■ j J d 3 *Li e* Xi ' ki <<5(k)d'(k / )5(xi) • • • <5(x n ))' c 

i= 1 J 


] [ / d 3 Xj e* Xi ki S n (k, xi, ■ ■ ■ x n ), 


(A.5) 


2=1 ' 


where 

S n (k;x i,---x n ) = (A.6) 

Note that the angle average is a linear operation and thus commutes with the Fourier trans¬ 
form; in other words, the £;-space angle average of the Fourier transform of a function is the 
Fourier transform of the x-space angle average of the same function. 
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Now consider the limit ki — > 0, which implies that Xi —> oo in the argument of S n . Then 
S n (k ) describes the modulation of the small-scale power spectrum P(k, 0 ) measured around 
x = 0 by n spherically symmetric large-scale modes (recall that k' ~ — k). This statement 
can be formalized by introducing an intermediate scale Rl such that 1/k <C Rl <C |xj| ~ 1 /ki 
and defining 5(k) — > 6r l ( k) to be the Fourier transform within a cubic volume of size Rl 
around x = 0. Then, J^ L (k) = <5(k) + 0(l/(kRL)), while the long-wavelength modes are 
constant over the same volume with corrections suppressed by kiRL • The corrections we 
expect in the end are thus of order ki/k. To lowest order in these corrections, S n {k ) can be 
written as 

lim : S n (k\ x\, ■ ■ ■ x n ) = ( ' 1 • • • [ (P(k, 0)<5 L (xi) • • • ^(x„))' . (A.7) 

ki —>-0 J 47r J 47r 

We can now insert the expression for the local power spectrum from Eq. (A.l), which imme¬ 
diately yields 

hin : S n (k;x i,---s n )= £ R m {k)P(k) j • f (C(0)5l( x i) • • • ^(x n ))'_ 0 . 

m =0 

(A. 8 ) 


Here, the subscript i — 0 indicates that only contractions between 0 and x ? ; are to be taken, 
since the l.h.s. is defined through the connected correlation function (all other contractions 
would contribute to the disconnected part of ((5(k)(5(k / )5(xi) • • • <5(x n ))). Since all density 
fields in the correlator in Eq. (A. 8 ) are linear, limiting to the contractions between 0 and x* 
then constrains rn = n. 4 We obtain 

1 TL n 

lim : S n (k\xi, ■ ■ ■ x n ) = — R n (k)P(k) n! TT / —— ^l(x*) 

0 n! - LJ - / 47 r 

i= 1 J 

'P' /» , 3 1 

= i 2 n (fc)P(A:) n J j^e^PUh). (A.9) 

Here, £l and Pl denote the linear matter correlation function and power spectrum, respec¬ 
tively. The angle average in the first line is trivial of course. Going back to Fourier space 
then immediately yields that for ki —> 0, 

S n (k-k 1 ,---k n ) = R n (k)P(k)f[P L (k i ) + o(^, , (A.10) 


where /cnl is the nonlinear scale. This can be reordered to yield Eq. (1.2), 


o /; \ 1 ■ S n {k]k\, ■ ■ ■ k n 

Rn(k) = hm — ;; -. 

P(k)P L (k!) • • • P L {k n ) 

This provides the connection between the response functions R n (k ) and the angle- 
averaged matter (n + 2)-point function Eq. (1.1) in a certain limit squeezed limit (since 
ki <C k). Note that no assumption about the magnitude of k has been made, i.e. this 
value can be fully nonlinear. In this paper, we accurately determine this fully nonlinear 
quantity using simulations. In the following subsections, we illustrate Eq. (1.2) at tree level 
in perturbation theory for the cases n = 1 (three-point function) and n = 2 (four-point 
function). 

4 If we use the definition of Eulerian response in Eq. (1.3) instead, we equivalently obtain a slightly more 
complicated relation where all m < n contribute. 
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A.l Tree-level result: n = 1 

At tree-level for n = 1 we obtain 


lim «Si(fc;lfei) tree — level 2 Km 

ki — >0 — ^0 


Eq. (1.2) then yields 


d 2 k 


Air 


- [E 2 (k, k JPlik) + F 2 {-k - ki, k^flk + ki|)] P t {k i) 


(A.ll) 


R l{ k) tvee J evel ^ hm J ^[F 2 (kM)W) + F2(-k-kiM)mk + ki\)] ■ (A. 12 ) 


where 


. 5 1 / k ki\ 2 2 

fMk.k,)(- + -A + -,. 


(A.13) 


where /x is the cosine of k and ki. The term n/2(k/ki) is problematic as we are sending 
k\ —> 0. Using that 


|k + ki| = k(l + qg + 0(q 2 )), q = y 


i^(|k + k 1 |)=^(*:) 


1 + ^T^ + °^ 


the sum of the two IR-divergent terms in Eq. (A. 12) becomes 

2 d In Pi(k) 


\ [4 P{k) + ~iktS> (|k+kil) } = l P k) 


p 


d In k 


- 1 


(A-14) 


+ 0{q). (A.15) 


As expected, the divergent pieces have canceled. We have dropped terms of order k\/k which 
are irrelevant in the limit we are interested in. We finally obtain 

r 1 


R\ (k) tree = evel 2 / ^ 

J-i 2 


TO 

1 

( 2 dlnP l . A . 4 2 

_ 7 

2 

dink + ) + 7^ J 


= 2 


10 ldlnPi 1 4' 

T ~~ 6 din A: ~~ 2 + 21 


47 1 d In P, 
21 3 d In k 


(A.16) 


This agrees with the linear prediction for Ri, obtained from substituting f\, e\, and g\ into 
Eq. (2.10). 


A.2 Tree-level result: n = 2 

At n = 2, we have 

d 2 k 2 


S n (k ; k\, k 2 ) = 

ki,k2 —>-0 


d 2 ki 


4-7T J 47T 
d 2 ki f d 2 k 2 


(5(k)<5(k / )5(k 1 )5(k 2 )) / c 


6 


F 3 (k,k 1 ,k 2 )P z (A : ) 


J 47T J 47T 
+ P 3 (-k- ki 2 ,ki,k 2 )Pz(|k + ki 2 |) 


+ 4 


P 2 (-k l5 k + k 1 )P 2 (k 2 , k + ki)P;(|k + kx|) 


+ P 2 (ki,k + k 2 )P 2 (-k 2 ,k + k 2 )P z (|k + k 2 |) ^(A; i)Pi(k 2 ), (A.17) 
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where F 3 is the symmetrized third-order perturbation theory kernel, and the unsymmetrized 
form can be found in [43]. Both F 2 and F 3 contain formally IR-divergent terms up to 
0 [(qi t 2 ) 2 } where q\;i = kip/k for k\p —> 0 which cancel in the end, so we should expand the 
power spectrum to 0 [(ki t 2 ) 2 ] to obtain the consistent result at order 0[{qi^) Q }- We have 


P l {\^ + \^ l \)=P l (k) 

Pzdk + kx + k.lHP^) 


Qiu 2i\ k dPi(k) qlul k 2 d 2 Pi(k) 3 


1 + ( gi/xi + —[1 - /xi ) )m) dk 


+ 


2 Pi(k) dk 2 
1 + ([<?iMi + Q 2 H 2 \ + J [qf( 1 ~ Mi) + <?!(! - Ml) 


+ <%?) 


+ 2gig 2 (Mi2 - M 1 M 2 )]) 


k dPi(k) 
Pl(k) d\nk 


+ \ (??Mi + q'Im! + 2qiq2di^2) d + 0[{ q i i2 ) 3 ] 


(A.18) 


where mi ,2 is the cosine of k and ki i2 , and M 12 is the cosine of ki and k 2 . The leading order 
terms are 


P 2 (fc) 


tree—level 


d 2 ki 

47T 


d 2 k 2 1 
4vr 147 


(628 + 324/i 2 + 112 mi 2 


— 28OM1M2M12 + 380m 2 + 656miM2 — 56m 2 Mi2) 


+ (—273mi + 147miM2Mi2 — 336m 2 — 483miM2 + 63m 2 Mi2) 


k dPi(k) 
Pi{k) dk 


+147miM2 


k 2 d 2 Pi{k) 
Pi{k) dk 2 


8420 100 k dPt(k) 1 k 2 d 2 Pi{k) 

1323 ~~ "63" Pi(k) dk + 9 P L (k) dk 2 ’ 


(A.19) 


which is in agreement with the linear prediction of Eq. (2.11) once the numbers for e*, /*, gi 
(i = 1,2) are inserted. 


B Analytical solution for S p and S a in Einstein-de Sitter 


As derived in [20], the modified cosmology described by a(t ) follows the evolution equation 
of a uniform density spherical perturbation in the fiducial cosmology. In this section, we 
take the fiducial cosmology to be Einstein-de Sitter (EdS), = 1, and consider a positive 
overdensity perturbation. In this case there is a well known parametric solution for the 
collapse [44]: 


d( 0 ) 

H 0 m 


1 Mr, 


2 tt rn ~ 1 


(1 — COS 0 ) 


n r 


2 (n m - 1 ) 3/2 


(9 — sin0). 


(B.l) 


Note that kl m — 1 = —£Ik > 0. We now want to derive the matching between Hq and 
with the fiducial values and the linearly extrapolated overdensity 5l o- The same matching 
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also works for 5lo < 0 in which case the sin, cos are to be replaced by sinh, cosh. First, the 
9 —> 0 limit should reduce to EdS, where 


a(t) = 



(B.2) 


so that a(to) = 1 and H(to) = a(to)/a(to) = Ho- One can further expand Eq. (B.l) around 
9 = 0. On the other hand, we know the linear relation between a(t ) and a(t), since at linear 
order 5 a = —5 p /3 = -5 L oa(t)/ 3: 


a(t) = a(t)[l + <5 a (t)] 





(B.3) 


where we have used that a(to) = 1 and 5lo is the linearly extrapolated overdensity at to- 
Matching at zeroth order in (Hot) 2 / 3 yields an identity. At linear order in 5lo > we obtain 


-= x</lo and 1 - (1 + 5h ) 2 = -Slo , 

l Z m O O 


(B.4) 


where 1 + 5h = Hq/Hq- Thus, we have determined all quantities in Eq. (B.l) in terms of 
<5lo- 

One might wonder what happens in the limit <5^0 —> 3/5, which implies K / Hq —> 1, 
—> oo and 5h —> — 1 and thus seemingly an ill-dehned cosmology. Let us investigate the 
solution Eq. (B.l) in this limit. First, we have 


Horn 


1 (i + M ~ 3 

2 (Q m - 1)3/2 


(6 — sin 9) 


^(1-(1 + 5 h ) 2 ) 3/2 (0-sin 9). 


(B.5) 


For Sh -> —1 and Q m oo, this solution becomes 


a(9) = -(1 — cos 9) 

Hot(9) = ^ (9 - sin 9). (B.6) 


Thus, the solution remains perfectly valid. It is simply the parametrization in terms of Hq 
and which breaks down. We also see why this happens: turn-around happens in Eq. (B.6) 
at 9 = 7T, Hotta = vr/2, and at a = 1- Thus, the Hubble constant goes to zero just at the point 
where we are trying to match the Friedmann equation, which of course assumes an expanding 
universe. For even larger overdensities, a(9) never reaches unity and thus no matching to 
Ho, is possible. Note that solutions of course exist, they are simply not captured by 
a parametrization of the form Eq. (B.l). In any case, such large nonperturbative values of 
the long-wavelength overdensity are not of practical interest for the application to separate 
universe simulations. 
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B.l Perturbative solution and nonlinear growth factor 

In this section we derive the series solution Eq. (2.7) for 5 a and 5 P in EdS. We write the 
parametric solution as 


d(9) = -e 1 (1 — cos0) 
to 4 

where we have used a(t ) = (t/to) 2 ^ 3 , to = 2/(377o), and defined 

~ “ 77( 1 5 


e = 


0, 


= 3&0 ■ 


Our goal is to obtain 
Thus, we need to solve 


a(t 0 ) = 1 + 5 a (to) • 


1 = t(0 0 ) <=> -e 3/2 = 9 0 - sin 9 0 


for 0q. We perform a series expansion, 

1 00 

@o - sm0 O = Q e 3 - m e 5 0 + ■■■ = £ M 0 2n+1 , 


71=1 


(B.7) 

(B.8) 

(B.9) 

(B.10) 

(B.ll) 


and solve Eq. (B.10) order by order, 
solution has to solve 


The leading solution is 0Q 1 ' = 2c 1 / 2 . The n-th order 


k =1 


q(n—fc+1) 
7 0 


2fc+l 


(B.12) 


Note that in order to trust the final expression at order 5™ 0 , this solution needs to be expanded 
to order n = m + 2. In the following, we choose m = 5. Solving this order by order, we 
obtain 


d 0 = 2c 1 / 2 


1 2 o 

1 + l5 £+ l75 e + 


1575 


3 4 3 4 

' + 67375 £ + 


(B. 13) 


We then insert this into a{6) and expand in e, replacing e with 5/3 5lo through Eq. (B.8). 
This yields 

OO 

1 + 5 a (to) = o,(9q) = 1 + 'y ] e n 5£g , (B.14) 

n= 1 


where the first few coefficients are 


1 1 23 1894 3293 

—; eo =-; e3 =-; 6 a = -; =- . 

3 21’ 1701 392931 1702701 


(B.15) 


We can generalize this to other times t by replacing 5lo with 5Loa(t), which is possible since 
EdS is scale free: 

OO 

5 a (t) = ^>2 e n [5 L0 a(t)] n . (B.16) 

n=1 


- 26 












Finally, we obtain the density at t through 


«5 p (t) = [l + 5 a (t)]- 3 -l, 


(B.17) 


and expanding in powers of 5lo■ This yields 


Sp(t) = ^2 fn[Swa{t)] n 
n= 1 

_ _ 17. 341 55805 _ 213662 

Ji- , J2 21 j /3 - 567 , U - 130977 > h- 729729’ 


(B.18) 


We have verified that Eq. (B.16) and Eq. (B.18) are accurate in ACDM as well when replacing 
a(t) with D(t)/D(to), where D(t) is the growth factor in the fiducial cosmology. 


C Linear growth in modified cosmology 


In this section we iteratively solve the growth factor D{t) in the modified cosmology, for a 
fiducial EdS background, to obtain a perturbative expansion in terms of 5lo [Eq. (2.14)]. 
The modified cosmology is described by 


H 2 (a) 


= m 




-—3 


T (1 tlm)^ 


-2 


= H 2 [a " 3 - 


ea 


-—21 


(C.l) 


where e is defined in Eq. (B.8) and we have used Cl m HQ = Q m HQ = Hq. Similarly, the mean 
background density in the curved universe is given by 

4t rG^(t) = ^L m Hla-\t) = ^H 2 a~ 3 (t). (C.2) 

On the other hand, in the background EdS cosmology we have 

H(a) = H 0 a~ 3/2 . (C.3) 


We want to derive the growth of density perturbations 5 S = p/p — 1, where the subscript 
s denotes that they are small-scale fluctuations, in order to distinguish from the density 
perturbations 5 P = p/p — 1 with respect to the EdS background. Note also that 5 S is defined 
with respect to the background density of the curved universe (i.e. we do not include the 
“reference density” contribution here). The growth equation is then given by 


$ a + 2H(t)5 a -4:irGp(t)6 8 = 0, 


(C.4) 


which becomes 

l + 2H 0 (a~ 3 - ea- 2 )^ 2 l - ^H 2 a~ 3 5 s = 0 

l s + 2i7 0 a" 3/2 (l - ea) l /H s - ^H 2 a~ 3 8 s = 0 . (C.5) 

Note that we have neglected the curvature contribution to the Poisson equation, which in¬ 
volves 3 K$> and the curved-space Laplacian. If K /~ 1, these terms only become relevant 
for small-scale modes that are of order the horizon. Since we are studying the subhorizon 
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evolution of small-scale modes, and moreover we always have K / H 2 } <C 1 for a flat fidu¬ 
cial cosmology, these terms are entirely negligible for our purposes. We now replace t with 
y = In a(t) as time coordinate, where a(t) is the scale factor in the EdS background. Dividing 
by H 2 , and inserting e = 5/35lq, we obtain 


dy 2 


d s + 


2A- 3 / 2 (l - jj<5 L0 aA a ) 


1/2 


3 

2 



= 0, 


(C.6) 


where we have defined A a (t) = 1 + 6 a (t). So far, everything is exact. We now perform a 
series expansion of Eq. (C.6) in 5lo- At zeroth order, A a = 1, and we obtain 


dy 2 


^ 0) + 


IAj(o) 

2 dy s 


-^ 0) = 0, 
2 


(C.7) 


which has growing and decaying modes of 5$ oc e y and e~ 3y ^ 2 . Thus, it is identical to the 
growth in the background EdS cosmology, as expected. In the following, we will drop the 
decaying mode following standard practice. Further, we will normalize 5^ to a(t) = e y at 
early times, and replace it with D to denote the small-scale growth factor. Thus, = 

a(t). 

The series expansion of A a is given in Eq. (B.16), which we can generalize to other 
times t by replacing Slo with 5Loa(t). We see that Eq. (C.6) can be expanded into a series 
in powers of 6lo a, leading to 


W> 6 + 


y! c m^L0 f 

_m= 0 


my 


dy 


E i ™ my 

a mOLO e 


_m =0 


D = 0, 


with coefficients c m , d m . Specifically, the coefficients are defined through 

/ r \ 1/2 o OO 

2 A - 3 / 2 1 - -5 L0 aA a - - = c m [S L oa(t)r 

' ' m= 0 

00 

9 a- 3 = 

m=0 

Correspondingly, we write the pure growing-mode solution as a series 


(C.8) 


(C.9) 


D(y) = J2gnS n L0 


n^ e (n+l) y 


n=0 


(C.10) 


with coefficients g n . Given our choice of normalization, we have go = 1. Thus, 

7 00 72 00 

-D(y) = + l)g n 5 n L0 e^ y - —D(y) = + 1 ?g n S n L ^ n+1)y . (C.ll) 

d n =0 d n=0 

Assuming we have a solution to order n — 1, the solution at order n then has to satisfy 


(n + l) 2 g n 5 n L0 e^ y + 22 g n - m 62o m e (n - m+ V y [( n-m+ 1 )c m 5^e my 

m —0 


dmSZbe 


m myi _ 


= 0. 


(C.12) 
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The time dependence e ny factors out, and we obtain a simple algebraic relation for g n in 
terms of {c m , d m , 5 m }”"l 1 0 : 

n 

(n + 1 ) 2 g n + ^2 9n—m [(n -m + 1 )c m - d m \ = 0 . (C.13) 

m =0 

For example, for n = 1 we have 


4^1 + [2c 0 - d 0 ] 5i + 5o [<h - di] = 0 , 

where Co = 1/2, do = 3/2, c\ = —2/3, di = 3/2. Thus, 

7 2 3 13 13 

2 y 3 2 6 y 21 


(C.14) 


(C. 15) 


This is the result for the leading 0(5lo) correction to the small-scale growth. Straightforward 
algebra yields the extension to higher order. We thus obtain 


D(t ) = a(t) 

where the first few coefficients are 


1 + ^>2 9n[hoa(t)}' 


n= 1 


5i = 


13 
21 ’ 


52 = 


71 

189’ 


53 = 


29609 

130977’ 


54 = 


691858 

5108103 


(C. 16) 


(C.17) 


In order to generalize from EdS to other cosmologies, we perform the usual replacement of 
a(t) —>■ D(t), where D(t ) is the growth factor in the fiducial cosmology normalized to a(t) 
during matter domination. Thus, we obtain 


D(t) = D(t ) \l + J2 


9n 


n=1 




m 

D(to)\ 


(C.18) 


D Transformation of power spectrum 

This section briefly derives the transformation of the power spectrum under a rescaling of 
spatial coordinates 

x = cx , (D-l) 

where c is a constant. This is necessary in order to calculate the full power spectrum response, 
since wavenumbers are defined in comoving coordinates and the modified scale factor a ^ a 
at fixed time. 

The correlation function of a scalar 5, which satisfies d(x) = d(x(x)), then transforms 
as (App. A in [45]) 


2c 


2c 


?(r) = (dUU(-h)> = (<M w)d(-^))=e(c- 1 f). 


(D.2) 


Note that for our purposes 5 transforms as scalar, since we take into account the change in 
the background reference density separately. Since 


tfc- 1 ?) = 


d 3 k 

( 27 t ) 3 


P(k) exp \ic ■ 


(D.3) 
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and P( k) is defined through 


1(f) 


d 3 k 

(27r) 3 


P(k) exp 


ik • r 


one immediately obtains 

d 3 kP(k)=[d 3 kP(k)\ k=cii . 


Since c is constant, this leads to 

P(k) = c 3 P(ck). 


(D.4) 


(D.5) 

(D.6) 
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